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Abstract 



Langevin simulation provides an effective way to study col- 
lisional effects in beams by reducing the six-dimensional 
Fokker-Planck equation to a group of stochastic ordinary 
differential equations. These resulting equations usually 
have multiplicative noise since the diffusion coefficients in 
these equations are functions of position and time. Con- 
ventional algorithms, e.g. Euler and Heun, give only first 
order convergence of moments in a finite time interval. In 
this paper, a stochastic leap-frog algorithm for the numeri- 
cal integration of Langevin stochastic differential equations 
with multiplicative noise is proposed and tested. The al- 
gorithm has a second-order convergence of moments in a 
finite time interval and requires the sampling of only one 
uniformly distributed random variable per time step. As 
an example, we apply the new algorithm to the study of a 
mechanical oscillator with multiplicative noise. 

1 INTRODUCTION 

Multiple Coulomb scattering of charged particles, also 
called intra-beam scattering, has important applications in 
accelerator operation. It causes a diffusion process of par- 
ticles and leads to an increase of beam size and emittance. 
This results in a fast decay of the quaUty of beam and re- 
duces the beam lifetime when the size of the beam is large 
enough to hit the aperture . 

An appropriate way to study the multiple Coulomb scat- 
tering is to solve the Fokker-Planck equations for the dis- 
tribution function in six-dimensional phase space. Never- 
theless, the Fokker-Planck equations are very expensive to 
solve numerically even for dynamical systems possessing 
only a very modest number of degrees of freedom. Trunca- 
tion schemes or closures have had some success in extract- 
ing the behavior of low-order moments, but the systematics 
of these approximations remains to be elucidated. On the 
other hand, the Fokker-Planck equations can be solved us- 
ing an equivalent Langevin simulation, which reduces the 
six-dimensional partial differential equations into a group 
of stochastic ordinary differential equations. Compared to 
the Fokker-Planck equation, stochastic differential equa- 
tions are not difficult to solve, and with the advent of mod- 
ern supercomputers, it is possible to run very large num- 
bers of realizations in order to compute low-order moments 
accurately. In general, the noise in these stochastic ordi- 
nary differential equations are multiplicative instead of ad- 
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ditive since the dynamic friction coefficient and diffusion 
coefficient in the Fokker-Planck equations depend on the 
spatial position. An effective numerical algorithm to inte- 
grate the stochastic differential equation with multiplicative 
noise will significantly improve the efficiency of large scale 
Langevin simulation. 

The stochastic leap-frog algorithms in the Langevin sim- 
ulation are given in Section II. Numerical tests of this algo- 
rithms is presented in Section III. A physical application 
of the algorithm to the multiplicative-noise mechanic os- 
cillator is given in Section IV. The conclusions are drawn 
in Section V. 

2 STOCHASTIC LEAP-FROG 
ALGORITHM 

In the Langevin simulation, the stochastic particle equa- 
tions of motion that follow from the Fokker-Planck equa- 



tion are (Cf. Ref. 



V, 

F 



DT{t), 



(1) 

(2) 



where F is the force including both the external force and 
the self-generated mean field space charge force, m is the 
mass of particle, is friction coefficient, D is the diffusion 
coefficient, and T{t) are Gaussian random variables with 
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In the case not too far from thermodynamic equilibrium, 
the friction coefficient is given as 
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and the diffusion coefficient D is D = vkT/m [^. Here, 
n(r) is the density of particle, r(r) is the temperature of of 
beam, Z is the charge number of particle, e is the charge of 
electron, A is the Coulomb logarithm, and k is the Boltz- 
mann constant. For the above case, noise terms enter only 
in the dynamical equations for the particle momenta. In 
Eqn. (|]) below, the indices are single-particle phase-space 
coordinate indices; the convention used here is that the odd 
indices correspond to momenta, and the even indices to the 
spatial coordinate. In the case of three dimensions, the dy- 
namical equations then take the general form: 
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±4 = ^4(2:3) 

2^5 = F5{xi,X2,X3,X4:,X5,Xg) + a55{x2,X4,X(i)S,5{t) 

xe = Fe(x5) (6) 

In the dynamical equations for the momenta, the first term 
on the right hand side is a systematic drift term which in- 
cludes the effects due to external forces and damping. The 
second term is stochastic in nature and describes a noise 
force which, in general, is a function of position. The noise 
^(t) is first assumed to be Gaussian and white as defined 
by Eqns. (3)-(Q). The stochastic leap-frog algorithm for 
Eqns. (H) is written as 
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The deterministic contribution Di {h) can be obtained us- 
ing the deterministic leap-frog algorithm. Here, the deter- 
ministic contribution Di{h) and the stochastic contribution 
Si{h) of the above recursion formula for one-step integra- 
tion are found to be 
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where Wi{h) is a series of random numbers with the mo- 
ments 

imh)) = imih))^) = ((Ti^.w)^) = (9) 

3 (10) 



Figure 1: Zero damping convergence test. {x'^{t)) at i = 6 
as a function of step size with white Gaussian noise. Solid 
fines represent quadratic fits to the data points (diamonds). 

where i? is a uniformly distributed random number on the 
interval (0,1). This trick significantly reduces the computa- 
tional cost in generating random numbers. 

3 NUMERICAL TESTS 

The above algorithm was tested on a one-dimensional 
stochastic harmonic oscillator with a simple form of the 
multipficative noise. The equations of motion were 



p = Fi{p,x) +aix)^{t) 
X ~ p 



(12) 
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This can not only be achieved by choosing true Gaussian 
random numbers, but also by using the sequence of random 
numbers following: 




R < 1/6 
l/6<i?<5/6 (11) 
5/6 < R 



where Fi{p,x) — —^p — rf'x and (j{x) — —ax. The 
stochastic leapfrog integrator for this case is given by 
Eqns. (||) (white noise) with the substitutions xi ~ p, 

X2 = X. 

As a first test, we computed {x'^) as a function of time- 
step size. To begin, we took the case of zero damping con- 
stant (7 = 0), where {x"^) can be determined analytically. 
The curve in Fig. |l] shows (a:^) at t = 6.0 as a function 
of time-step size with white Gaussian noise. Here, the pa- 
rameters f] and a are set to 1.0 and 0.1. The analytically 
determined value of (x^) at t = 6.0 is 2.095222. The 
quadratic convergence of the stochastic leap-frog algorithm 
is clearly seen in the numerical results. We also verified 
that the quadratic convergence is present for nonzero damp- 
ing (7 = 0.1). Alt — 12.0, and with all other parameters 
as above, the convergence of {x'^) as a function of time step 
is shown by the curve in Fig. ^ As a comparison against 
the conventional Heun's algorithm [^, we computed (x^) 
as a function of t using 100, 000 numerical realizations for 
a particle starting from (0.0, 1.5) in the {x,p) phase space. 
The results along with the analytical solution and a numer- 
ical solution using Heun's algorithm are given in Fig. ||. 
Parameters used were h = 0.1, 77 = 1.0, and a = 0.1. The 
advantage in accuracy of the stochastic leap-frog algorithm 
over Heun's algorithm is clearly displayed, both in terms 
of error amplitude and lack of a systematic drift. 





Figure 2: Finite damping (7 = 0.1) convergence test. 
{x^it)) at t = 12 as a function of step size with white 
Gaussian noise. Solid lines represent quadratic fits to the 
data points (diamonds). 
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Figure 3: Comparing stochastic leap-frog and the Heun al- 
gorithm: (x^ {t)) as a function of t. Errors are given relative 
to the exact solution. 



4 APPLICATION 

In this section, we apply our algorithm to studying the ap- 
proach to thermal equilibrium of an oscillator with multi- 
plicative noise. The governing equations are: 
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where the diffusion coefficients D = XkT, A is the cou- 
pling constant, and luq is the oscillator angular frequency 
without damping. In Fig. ^ we display the time evolu- 
tion of the average energy with multiplicative noise from 
the simulations and the approximate analytical calcula- 
tions [^. The analytic approximation resulting from the 
application of the energy-envelope method is seen to be in 
reasonable agreement with the numerical simulations for 
kT = 4.5. The slightly higher equilibrium rate from the 



Figure 4: Temporal evolution of the scaled average energy 
{E{t)) with multiplicative noise from numerical simulation 
and analytical approximation. 



analytical calculation is due to the truncation in the energy 
envelope equation using the (£'^(t)) « 2{E{t))'^ relation 
which yields an upper bound on the rate of equilibration of 
the average energy . 

5 CONCLUSIONS 

We have presented a stochastic leap-frog algorithm for 
Langevin simulation with multiplicative noise. This 
method has the advantages of retaining the symplectic 
property in the deterministic limit, ease of implementa- 
tion, and second-order convergence of moments for mul- 
tiplicative noise. Sampling a uniform distribution instead 
of a Gaussian distribution helps to significantly reduce the 
computational cost. A comparison with the conventional 
Heun's algorithm highlights the gain in accuracy due to the 
new method. Finally, we have applied the stochastic leap- 
frog algorithm to a nonlinear mechanic-oscillator system to 
investigate the the nature of the relaxation process. 
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